In this analysis, explores spatial and raster data for the Mediterranean region, with a focus on Greece. The objective is to visualize and analyze elevation data within the administrative boundaries of Greece and associate elevation information with key cities.
Load Libraries:
library(sf) # classes and functions for vector data
library(terra) # classes and functions for raster data
library(spData) # load geographic data
library(spDataLarge) # load larger geographic data
library(dplyr)
library(mapsf)
library(tmap)
library(leaflet) # for interactive maps
library(ggplot2)
Load Raster Data :Loading and plotting a raster file containing elevation data for the Mediterranean region.
raster_file = rast("30n000e_20101117_gmted_mea300.tif")
class(raster_file)
## [1] "SpatRaster"
## attr(,"package")
## [1] "terra"
plot(raster_file, col = terrain.colors(6))
There are some faint distinctions of Greece to the right, Italy in the center, and above it, the Alps visible towards the upper left.
Read Shapefiles : Reading three shapefiles: one for the administrative boundaries of Greece (GRC_adm2.shp), one for their capitals (poleis.shp), and one for general locations (places.shp).
nomoi <- st_read("GRC_ADM2/GRC_adm2.shp")
## Reading layer `GRC_ADM2' from data source
## `C:\Users\Stella\OneDrive\Υπολογιστής\[02] Διερευνητική Ανάλυση και Οπτικοποίηση Δεδομένων\hw4\data μέρος β\GRC_ADM2\GRC_ADM2.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 77 features and 7 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 19.3726 ymin: 34.8009 xmax: 29.60684 ymax: 41.74889
## Geodetic CRS: WGS 84
poleis <- st_read("poleis/poleis.shp")
## Reading layer `poleis' from data source
## `C:\Users\Stella\OneDrive\Υπολογιστής\[02] Διερευνητική Ανάλυση και Οπτικοποίηση Δεδομένων\hw4\data μέρος β\poleis\poleis.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 51 features and 2 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 149581.4 ymin: 3895170 xmax: 878802.9 ymax: 4555185
## Projected CRS: GGRS87 / Greek Grid
places <- st_read("places/places.shp")
## Reading layer `places' from data source
## `C:\Users\Stella\OneDrive\Υπολογιστής\[02] Διερευνητική Ανάλυση και Οπτικοποίηση Δεδομένων\hw4\data μέρος β\places\places.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 4326 features and 4 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 19.39397 ymin: 34.84761 xmax: 27.2865 ymax: 41.72253
## Geodetic CRS: WGS 84
An important observation is that the nomoi file comprises 77 administrative regions, whereas the capitals consist of 51. This disparity will be revisited later. Furthermore, it’s noteworthy that the geometry of nomoi is multipolygon using longitude and latitude, while the capitals and places have point geometries and are projections in Greek Grid. This aspect is addressed further down in the analysis to ensure a cohesive representation when plotting them on the same map.
Crop Raster to Greece’s Bounding Box:
Extracting the bounding box (the rectangular area that contains all the points) of the administrative boundaries of Greece.
# Get the bounding box of Greece
bbox_greece <- st_bbox(nomoi)
Using this bounding box to crop the elevation data raster to include only the area of Greece.
# Crop raster to the bounding box of Greece
raster_greece <- crop(raster_file, bbox_greece)
Mask Raster with Greece’s Geometry:
Using the actual geographical shape of Greece to mask the elevation data raster. This means keeping only the elevation values that fall within the borders of Greece.
# Mask raster using the geometry of Greece
raster_greece <- mask(raster_greece, nomoi)
Plot the resulting raster: Finally, plotting the newly created raster file, which now contains elevation data only for Greece.
plot(raster_greece, col = terrain.colors(6))
The towering peak of Mount Olympus, the highest point in Greece, prominently marks the map with its distinctive white hue. Other smaller white points, though present, are less conspicuous due to their limited area size.
City Elevation Analysis: To accurately depict city locations on the regional map, it is imperative to transform the city coordinate reference system to match the longitude-latitude format used for the regions:
poleis_latlon <- st_transform(poleis, crs(raster_greece))
Then calculating elevation values for each city(regional capitals) within Greece:
# Extract elevation values from the raster map at the locations of cities
elev_values <- extract(raster_greece, poleis_latlon)
# Add elevation values in poleis_latlon
poleis_latlon$elevation <- elev_values
Visualization of Administrative Boundaries and Cities with Elevation Information: Finally, plotting the administrative boundaries of Greece and overlay capitals with elevation information, represented by red points whose size is proportional to the elevation:
# Extract the elevation and normalize it
cex <-poleis_latlon$elevation[2] / (max(poleis_latlon$elevation)/2)
# Plot the raster map of Greece
plot(raster_greece, col = terrain.colors(6), reset = FALSE)
# Plot the boundaries of the regions ("nomoi") on top of the raster map
plot(nomoi["geometry"],add = TRUE,axes = TRUE)
# Plot the red points representing the regional capitals on top of the map
# Use the normalized elevation value to adjust the size of the points
plot(poleis_latlon[0], add = TRUE, col = "red", pch = 16, cex = cex[,1])
(Some points may be challenging to observe due to their smaller size in comparison to the larger ones)
The majority of higher-elevated cities are concentrated in the northern regions of the country. Among the capitals, Karpenisi boasts the highest elevation. This finding aligns with the known geography of central Greece. Nestled amidst mountainous landscapes, as illustrated in the elevation map of Greece, Karpenisi is situated in proximity to the Pindos mountain range. Notably, the capitals situated near the Pindos mountain range exhibit higher elevations, whereas those situated near the sea generally attain lower elevations. Reinforcing this, the capital with the lowest elevation is Leykada.
Moreover, it is probable that regions situated near the sea would instinctively position their capitals near ports, whereas inland areas, lacking direct access to the sea, may have historically chosen elevated locations for enhanced protection. This strategic decision-making traces back to ancient times, contributing to the development of these areas that were eventually established as capitals in the modern era.
Data Preparation:
# Extract elevation values for each region
elev_values <- extract(raster_greece, nomoi)
# Factorize the ID column to perform a group by on it
elev_values$ID <- factor(elev_values$ID)
# Rename the elevation column
elev_values <- elev_values %>%
rename_at(vars(2), ~ "elevation")
# Group by the "ID" column and Summarize the grouped data (calculate mean and sd)
elev_values <- elev_values %>% group_by(ID) %>% summarize(mean_value = mean(elevation),sd_value = sd(elevation))
# Add mean and sd to the original nomoi dataset
nomoi$mean_elevation <- elev_values[[2]]
nomoi$sd_elev <- elev_values[[3]]
The dataset nomoi now includes two new columns, mean_elevation and sd_elev representing the mean and standard deviation of elevation for each region.
Visualizing Mean Elevation
Visualizing the geographic distribution of mean elevation across
administrative regions in Greece:
mf_theme("default")
mf_map(
x = nomoi,
var = "mean_elevation",
type = "choro",
breaks = "geom",
nbreaks = 5,
pal = terrain.colors(8, alpha=1, rev = FALSE),
border = "white",
lwd = 0.5,
leg_pos = "topright",
leg_title = "Mean Elevation",
)
The choropleth map provides a visual representation of the mean elevation across different administrative regions in Greece. Regions with higher mean elevations are depicted with distinct colors, allowing for a quick and intuitive understanding of the geographic patterns.
Comments:
As anticipated, the highest mean elevation is notably found in the areas
surrounding the Pindos mountain range. As one moves away from the
central mainland towards the coast, there is a discernible decrease in
the mean elevation.
Two intriguing exceptions to this trend are the islands of Crete and Ikaria, that stand out for their remarkable mountainous terrain, defying the typical expectations for an island.
Visualizing Standard Deviation of Elevation
Similarly, a choropleth map to visualize the standard deviation of
elevation across administrative regions:
mf_theme("default")
mf_map(
x = nomoi,
var = "sd_elev",
type = "choro",
breaks = "geom",
nbreaks = 5,
pal = terrain.colors(8, alpha=1, rev = FALSE),
border = "white",
lwd = 0.5,
leg_pos = "topright",
leg_title = "Standar Deviation of Elevation"
)
Comments:
Once again, the Pindos mountain range stands out as the focal point
where elevated standard deviation values are consistently noted. This
observation implies a notably diverse topography within these
regions.
In fact, the majority of Greece showcases a varied landscape, featuring a mix of elevated and lower terrains across its expanse. This characteristic is emblematic of a mountainous terrain enveloped by the sea, with coastlines in close proximity, contributing to the rich topographical tapestry of the country.
Spatially associating each region with its corresponding capital, determining the region to which each capital belongs. Then, computing the absolute difference between the mean elevation of the region and the elevation of its capital.
# Ensure the validity of administrative region geometries
nomoi <- st_make_valid(nomoi)
# Spatially join cities with administrative regions based on the "within" relationship
nomoi_with_capitals <- st_join(poleis_latlon,nomoi, join = st_within)
# Calculate absolute difference in elevation between cities and their associated administrative regions
nomoi_with_capitals$abs_diff_elevation <- unlist(abs(nomoi_with_capitals$elevation[2] - nomoi_with_capitals$mean_elevation))
Data preprocessing: As discovered previously the regions in the data are more than the capitals. Thus, filtering the dataframe nomoi to keep only the regions with capitals:
# Find the intersection of feature_id between nomoi and nomoi_with_capitals
inter <- intersect(nomoi$feature_id,nomoi_with_capitals$feature_id)
# Filter nomoi based on the common feature_id values obtained from the intersection
filtered_nomoi <- nomoi %>%
filter(feature_id %in% inter)
# Add a new column 'abs_diff_elevation' to filtered_nomoi and populate it with corresponding values from nomoi_with_capitals
filtered_nomoi$abs_diff_elevation <- nomoi_with_capitals$abs_diff_elevation[match(filtered_nomoi$feature_id, nomoi_with_capitals$feature_id)]
The extra regions are just further sub divisions of larger regions (e.g. Attica -> West,North,South,East Attica) or errors:
diff <- nomoi %>%
filter(feature_id %in% setdiff(nomoi$feature_id,nomoi_with_capitals$feature_id))
diff$NAME
## [1] "Regional Unit North Athens" "Regional Unit of West Attica"
## [3] "Mykonos Regional Unit" "Karpathos Regional Unit"
## [5] "???????????? ??????? ???? - ??????" "???????????? ??????? ????????"
## [7] "Regional Unit of Lemnos" "Regional Unit of Ikaria"
## [9] "Regional Unit of Piraeus" "Regional Unit of West Athens"
## [11] "???????????? ??????? ??????" "???????????? ??????? ??"
## [13] "???????????? ??????? ???????" "Thira Regional Unit"
## [15] "Regional Unit of Sporades" "Paros Regional Unit"
## [17] "Milos Regional Unit" "Milos Regional Unit"
## [19] "Regional Unit of East Attica" "Naxos Regional Unit"
## [21] "Thasos Regional Unit" "Ithaca Regional Unit"
## [23] "???????????? ??????? ?????" "Regional Unit of South Athens"
## [25] "???????????? ??????? ?????" "Regional Unit of Islands"
Visualizing Absolute Difference of Elevation: A choropleth map to visualize the absolute difference of elevation between the regions and their capitals:
mf_theme("default")
# plot population density
mf_map(
x = filtered_nomoi,
var = "abs_diff_elevation",
type = "choro",
breaks = "geom",
nbreaks = 5,
pal = terrain.colors(8, alpha=1, rev = FALSE),
border = "white",
lwd = 0.5,
leg_pos = "topright",
leg_title = "Absolute Difference"
)
The substantial variance between the mean elevation of regions and the elevations of their respective capitals is prevalent across many areas. While most capitals are situated near sea level and boast relatively low altitudes, the presence of towering mountains in these regions significantly raises the overall mean elevation. This phenomenon results in a pronounced disparity. Even capitals positioned away from coastal areas exhibit a substantial deviation from the highest points within the region. With the exception only of a handful of capitals which align with the general altitudes prevalent in their respective regions.
The top 10 regions based on the mean elevation:
nomoi[order(nomoi$mean_elevation, decreasing = TRUE), "NAME"][1:10,]$NAME
## [1] "Kastoria Regional Unit" "Regional Unit of Evrytania"
## [3] "Florina Regional Unit" "Ioannina Regional Unit"
## [5] "Grevena Regional Unit" "Regional Unit of Phocis"
## [7] "Regional Unit of Kozani" "Arcadia Regional Unit"
## [9] "Regional Unit of Trikala" "Drama Regional Unit"
The top 10 regions based on the standard deviation of elevation:
nomoi[order(nomoi$sd_elev, decreasing = TRUE), "NAME"][1:10,]$NAME
## [1] "Imathia Regional Unit" "Pieria Regional Unit"
## [3] "Regional Unit of Trikala" "Pella Regional Unit"
## [5] "Regional Unit of Phocis" "Chania Regional Unit"
## [7] "Arta Regional Unit" "Achaea Regional Unit"
## [9] "Karditsa Regional Unit" "Ioannina Regional Unit"
Preprocessing:
Converting the reference system of places to longitude-latitude:
places_latlon <- st_transform(places, crs(raster_greece))
Calculating elevation values for each place within Greece:
elev_values <- extract(raster_greece, places_latlon)
colnames(elev_values) <- c("ID","value")
# Add elevation values in places_latlon
places_latlon$elevation <- elev_values
Upon examination, it becomes evident that certain data points exhibit discrepancies; they either represent errors, given that they do not correspond to locations within Greece, or they are absent from the regional data used to derive the elevation raster for Greece. Consequently, the elevation values for these data points are recorded as NA:
places_latlon[is.na((places_latlon[["elevation"]][,2])),]
## Simple feature collection with 100 features and 5 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 19.49628 ymin: 35.61053 xmax: 27.04166 ymax: 41.72253
## Geodetic CRS: WGS 84
## First 10 features:
## osm_id name type population
## 89 61886594 Astakida island 0
## 107 65035962 Ψυττάλεια (Psyttalia) island 0
## 119 83942427 Pelouzo island 0
## 120 83942679 Pontikonisia island 0
## 124 83947121 Άγιος Ιωάννης (Agios Ioannis) island 0
## 478 258063721 Gramvousa Island island 0
## 669 273814587 Генералово village 0
## 675 273874873 Горна Арда village 0
## 676 273877347 Капитан Андреево village 0
## 785 282497732 Marathonisi island 0
## geometry elevation.ID elevation.value
## 89 POINT (26.82026 35.88547) 89 NA
## 107 POINT (23.58804 37.94123) 107 NA
## 119 POINT (20.94805 37.70375) 119 NA
## 120 POINT (20.8625 37.68604) 120 NA
## 124 POINT (20.62135 37.8204) 124 NA
## 478 POINT (23.57941 35.61053) 478 NA
## 669 POINT (26.27687 41.72253) 669 NA
## 675 POINT (24.60772 41.44903) 675 NA
## 676 POINT (26.32017 41.72034) 676 NA
## 785 POINT (20.86908 37.68582) 785 NA
Among these locations, the ones which are located within Greece exclusively comprise small islands, as discerned from their zero populations, unlikely to surpass 1500 meters. Consequently, these islands can be confidently omitted from consideration.
# Filter places_latlon to include only locations with elevation higher than 1500 meters
places_higher_than_1500 <- na.omit(places_latlon[places_latlon[["elevation"]][,2] > 1500,])
# Keep only inhabited places (population > 0) among those with elevation higher than 1500
inhabited <- subset(places_higher_than_1500, population > 0)
# Keep only uninhabited places (population == 0) among those with elevation higher than 1500
uninhabited <- subset(places_higher_than_1500, population == 0)
inhabited
## Simple feature collection with 1 feature and 5 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 21.01633 ymin: 40.10259 xmax: 21.01633 ymax: 40.10259
## Geodetic CRS: WGS 84
## osm_id name type population geometry
## 1489 296596658 Σαμαρίνα hamlet 701 POINT (21.01633 40.10259)
## elevation.ID elevation.value
## 1489 1489 1533
uninhabited
## Simple feature collection with 6 features and 5 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 21.10604 ymin: 39.576 xmax: 21.31684 ymax: 40.84674
## Geodetic CRS: WGS 84
## osm_id name type population geometry
## 1174 286918633 Αθαμανία region 0 POINT (21.18557 39.576)
## 1177 287220488 Oros Triklario region 0 POINT (21.10604 40.67903)
## 1190 287224403 Oros Varnoundas region 0 POINT (21.23711 40.84674)
## 1191 287224488 Vitsi peak 0 POINT (21.31684 40.72146)
## 1192 287224694 Verno (Vitsi) region 0 POINT (21.31392 40.75021)
## 2998 344746654 Λάκμος (Περιστέρι) region 0 POINT (21.114 39.67118)
## elevation.ID elevation.value
## 1174 1174 1759
## 1177 1177 1606
## 1190 1190 1969
## 1191 1191 1542
## 1192 1192 1581
## 2998 2998 1964
tmap_mode("view")
names(raster_greece) <- "elevation"
tm_shape(raster_greece) + tm_raster("elevation", palette = terrain.colors(6)) + tm_shape(nomoi) + tm_borders() + tm_shape(uninhabited) + tm_dots(col = "red", size = 0.1) +tm_text("name", just = "top") + tm_shape(inhabited) + tm_dots(col = "green", size = 0.1) + tm_text("name", just = "top")